Multi-terminal differential protection setting method for distributed generator t-connected distribution network

ABSTRACT

The present invention proposes a multi-terminal differential protection setting method for a distributed generator T-connected distribution network. The present invention can accurately identify internal faults and external faults of the protected area in the high-proportion distributed generator T-connected distribution network, with sufficient reliability and sensitivity.

TECHNICAL FIELD

The present invention relates to the technical field of electrical technologies, and particularly relates to a multi-terminal differential protection setting method for a distributed generator T-connected distribution network.

BACKGROUND

Building a new power system based on clean and low-carbon energy is a key measure to implement the two-carbon strategy. Large-scale multi-point connection of a distributed generator into a distribution network is one of the most important forms of new energy interconnection. In the vast rural distribution networks, inverter interfaced distributed generators (IIDGs) such as photovoltaic type and energy storage type are basically connected into the distribution network by means of T-connection. How to realize the differential protection after multi-point T-connection of high-proportion IIDGs into the distribution network has drawn much attention in the industry. For a distribution network containing T-connected IIDGs, two-terminal differential protection usually includes the current of T-connected lines into the protection action threshold, and the characteristics of the protection action are affected by the capacity of the T-connected IIDGs. When the capacity of the connected IIDGs is large, refused operation of fault protection will occur in the area. With the development of communication technologies and the wide application of primary and secondary fusion intelligent switches, multi-terminal differential protection has become a feasible and even preferable solution. In terms of the principle and scheme of multi-terminal differential protection for a distribution network applicable to multi-point T-connected lines, some scholars have proposed a multi-point T-connected line differential protection scheme based on master-slave topology, which can accurately identify internal and external faults of the area to a certain extent, but does not fully consider the influence of fault characteristics of the IIDGs on the amplitudes and phases of differential current and braking current. Some scholars have proposed a differential protection scheme based on the amplitude-phase relationship of multi-terminal current, which considers the influence of fault characteristics of the IIDGs on differential current, but does not fully consider the influence of the data synchronization error on differential current. Some scholars have studied the possibility of refused operation and maloperation of multi-terminal differential protection caused by the data synchronization error based on the engineering example of T-connected lines in the distribution network, which proves that the data synchronization error will have a certain effect on the performance of multi-terminal differential protection.

In the prior art, internal faults and external faults of the protected area cannot be accurately identified in the high-proportion distributed generator T-connected distribution network, the calculated result of the differential current of superimposed multi-terminal current phasor errors is not accurate, and the speed for calculating the overall maximum error of differential current of superimposed multi-terminal current phasor errors is low.

SUMMARY

To solve the above problem, the purpose of the present invention is to propose a multi-terminal differential protection setting method for a distributed generator T-connected distribution network, which can accurately identify internal faults and external faults of the protected area in the high-proportion distributed generator T-connected distribution network, with sufficient reliability and sensitivity, and can quickly and accurately calculate the overall maximum error of differential current of superimposed multi-terminal current phasor errors.

To achieve the above technical purpose, the present invention provides a multi-terminal differential protection setting method for a distributed generator T-connected distribution network, comprising the following steps:

S1: Connecting a distributed generator into a distribution network by means of T-connection, and making the distributed generator equivalent to a positive sequence current source controlled by grid interconnection voltage;

S2: Acquiring multi-terminal current according to the positive sequence current source controlled by grid interconnection voltage;

S3: Calculating the overall maximum error of superimposed multi-terminal current errors according to the multi-terminal current;

S4: Establishing non-convex constraint space according to the overall maximum error of superimposed multi-terminal current errors;

S5: Reducing and discretizing the non-convex constraint space;

S6: Obtaining the optimization problem of the non-convex constraint space based on the reducing and discretizing results of the non-convex constraint space;

S7: Solving the optimization problem of the non-convex constraint space;

S8: Analyzing and verifying the solution result of the optimization problem of the non-convex constraint space.

Preferably, in step S1, an equivalent model for making the distributed generator equivalent to a positive sequence current source controlled by grid interconnection voltage is:

$\begin{matrix} {{U_{{pcc}(1)}^{*} = {U_{{pcc}(1)}/U_{{{pcc}(1)}.N}}}{I_{{DGm}.N} = {\sqrt{2}P_{{ref}.N}/3U_{{{pcc}(1)}.N}}}{i_{q(1)} = \left\{ {{\begin{matrix} {0,} & {U_{{pcc}(1)}^{*} > 0.9} \\ {{{- 2}\left( {1 - U_{{pcc}(1)}^{*}} \right)I_{{DGm}.N}},} & {0.4 \leq U_{{pcc}(1)}^{*} \leq 0.9} \\ {{{- 1.2}I_{{DGm}.N}},} & {U_{{pcc}(1)}^{*} < 0.4} \end{matrix}i_{d(1)}} = {{{\min\left( {\frac{\sqrt{2}P_{ref}}{3U_{{pcc}(1)}},\sqrt{\left( {1.2I_{{DGm}.N}} \right)^{2} - i_{q(1)}^{2}}} \right)}i_{DG}} = \left\{ \begin{matrix} {{\frac{\sqrt{i_{d(1)}^{2} + i_{q(1)}^{2}}}{\sqrt{2}}{\angle\left( {{{arc}\tan\frac{i_{q(1)}}{i_{d(1)}}} + \varphi} \right)}},} & {i_{d(1)} \neq 0} \\ {{\frac{\sqrt{i_{d(1)}^{2} + i_{q(1)}^{2}}}{\sqrt{2}}{\angle\left( {{{- 90}{^\circ}} + \varphi} \right)}},} & {i_{d(1)} = 0} \end{matrix} \right.}} \right.}} & (1) \end{matrix}$

wherein U_(pcc(1)) is an effective value of grid interconnection positive sequence phase voltage, U_(pcc(1).N) is an effective value of grid interconnection positive sequence phase voltage in the rated running state, I_(DGm.N) is an amplitude of rated current, P_(ref.N) is a rated active reference value, P_(ref) is an active reference value before failure, İ_(DG) is output current, {dot over (ι)}_(d(1)) and {dot over (ι)}_(q(1)) are a d-axis component and a q-axis component of the positive sequence component of the output current, and φ is the phase of {dot over (U)}_(pcc(1)).

Preferably, step S3 of calculating the overall maximum error of superimposed multi-terminal current errors according to the multi-terminal current comprises:

Expressing the power frequency current phasor involved in a group of n(n≥3)-terminal current differential protection by a set P:

P={İ ₁ ,İ ₂ , . . . ,İ _(n)}  (2)

The action criterion of the multi-terminal current differential protection is:

I _(op) ≥I _(set)  (3)

wherein I_(op) is differential current; and I_(set) is a protection setting value;

Expressing the relationship between the measured value and the actual value of each power frequency current phasor as follows:

İ _({dot over (ι)}.ms) =İ _({dot over (ι)}.tr) +Δİ _({dot over (ι)})  (4)

i=1, 2, . . . , n

wherein İ_({dot over (ι)}.ms) represents the measured value of the power frequency current phasor, İ_({dot over (ι)}.tr) represents the actual value of the power frequency current phasor, and Δİ_({dot over (ι)}) represents the error between the actual value and the measured value;

At this moment, the differential current is:

$\begin{matrix} \begin{matrix} {I_{op} = {❘{{\overset{.}{I}}_{1.{ms}} + {\overset{.}{I}}_{2.{ms}} + \ldots + {\overset{.}{I}}_{n.{ms}}}❘}} \\ {= {❘{{\overset{.}{I}}_{1.{tr}} + {\Delta{\overset{.}{I}}_{1}} + {\overset{.}{I}}_{2.{tr}} + {\Delta{\overset{.}{I}}_{2}} + \ldots + {\overset{.}{I}}_{n.{tr}} + {\Delta{\overset{.}{I}}_{n}}}❘}} \end{matrix} & (5) \end{matrix}$

According to the Kirchhoff s Current Law, meeting the formula below during normal operation and during failure outside a protected area:

İ _(1.tr) +İ _(2.tr) + . . . +İ _(n.tr)=0  (6)

According to formula (5) and formula (6), the differential current during normal operation and during failure outside a protected area is:

I _(op) =|Δİ ₁ +Δİ ₂ + . . . +Δİ _(n)|  (7)

Denoting the overall maximum error of superimposed multi-terminal current errors as ΔI_(Σ):

ΔI _(Σ) =|Δİ ₁ +Δİ ₂ + . . . +Δİ _(n)|  (8)

Preferably, step S4 of establishing non-convex constraint space according to the overall maximum error of superimposed multi-terminal current errors comprises:

According to Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n), establishing a sector shadow area, wherein the phasors of Δİ*₁,Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n) end on the boundaries of the respective sector shadow area;

If the position relationship between Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n) in a complex plane is that Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n) are at an angle of θ*_(n);

wherein θ*_(n) is the phase difference between phasors Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n), and θ*_(n)∈[0, 2π);

According to the Law of Cosines, ΔI_(Σ.max) ² is:

$\begin{matrix} {{\Delta I_{\sum{.\max}}^{2}} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}}}❘}^{2} + {❘{\Delta{\overset{.}{I}}_{n}^{*}}❘}^{2} - {{2 \cdot {❘{\Delta{\overset{.}{I}}_{n}^{*}}❘} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}}}❘}}{\cos\left( {\pi - \theta_{n}^{*}} \right)}}}}} & (14) \end{matrix}$

Denoting Δİ*_(n) having the phase difference of θ*_(n) from Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) as Δİ_(n) ^(θ) ^(n) ^(*), and denoting Δİ_(n) ^(θ) ^(n) ^(*) ending on the boundaries of the shadow area as Δİ_(n) ^(θ) ^(n) ^(*,max);

Since Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ_(n) ^(θ) ^(n) ^(*) are independent of each other, the derivative of ΔI_(Σ) ² with respect to |Δİ_(n) ^(θ) ^(n) ^(*)| is:

$\begin{matrix} {\frac{d\left( {\Delta I_{\sum}^{2}} \right)}{d\left( {❘{\Delta{\overset{.}{I}}_{n}^{\theta_{n}^{*}}}❘} \right)} = {{2{❘{\Delta{\overset{.}{I}}_{n}^{\theta_{n}^{*}}}❘}} - {2 \cdot {\cos\left( {\pi - \theta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}}}❘}}}} & (15) \end{matrix}$

The second derivative of ΔI_(Σ) ² with respect to |Δİ_(n) ^(θ) ^(n) ^(*)| is:

$\begin{matrix} {\frac{d^{2}\left( {\Delta I_{\sum}^{2}} \right)}{{d\left( {❘{\Delta{\overset{.}{I}}_{n}^{\theta_{n}^{*}}}❘} \right)}^{2}} = {2 > 0}} & (16) \end{matrix}$

According to formula (14) to formula (16), since the quadratic term coefficients of the second derivative and the first derivative are both larger than 0, as |Δİ_(n) ^(θ) ^(n) ^(*)| changes from 0 to |Δİ_(n) ^(θ) ^(n) ^(*, max)|, ΔI_(Σ) has at most one extreme point which is a minimum point;

ΔI_(Σ) is maximum when Δİ_(n) ^(θ) ^(n) ^(*)=0 or Δİ_(n) ^(θ) ^(n) ^(*)=Δİ_(n) ^(θ) ^(n) ^(*, max);

At least one Δİ_(n) ^(ξ) inevitably exists in the sector shadow area so that:

|Δİ* ₁ +Δİ* ₂ + . . . +Δİ _(n-1) +Δİ _(n) ^(ξ) |>|Δİ* ₁ +Δİ* ₂ + . . . +Δİ* _(n-1)+0|  (17)

When the overall error of superimposed multi-terminal current errors is maximum, Δİ*_(n)=Δİ_(n) ^(θ) ^(n) ^(*, max);

Numbering four boundaries of the sector shadow area corresponding to Δİ_(n),

The boundaries 1 and 2 are straight boundaries, and the boundaries 3 and 4 are arc boundaries;

wherein Δİ_(n) ending on the boundary 1 can be decomposed into:

Δİ _(n) =Δİ* _(n,x) +Δİ _(n,y)  (18)

wherein Δİ*_(n,x) is the projection of Δİ_(n) on the vertical line of the boundary 1, and Δİ_(n,y) is the projection of Δİ_(n) on the boundary 1;

When Δİ*_(n,y) ends near the A₁ side, ΔI_(Σ.max) ² is:

$\begin{matrix} {{\Delta I_{\sum{.\max}}^{2}} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}} + {\Delta{\overset{.}{I}}_{n,y}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}^{2} + {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘}^{2} - {{2 \cdot {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}{\cos\left( {\frac{\pi}{2} - \delta_{n}^{*}} \right)}}}}}} & (19) \end{matrix}$

wherein δ*_(n) is the phase difference between Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n,x), and δ*_(n) ∈[0, 2π);

Since Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n,x), Δİ_(n,y) are independent of each other, the derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is:

$\begin{matrix} {\frac{d\left( {\Delta I_{\sum}^{2}} \right)}{d\left( {❘{\Delta{\overset{.}{I}}_{n,y}}❘} \right)} = {{2{❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘}} - {2{{\cos\left( {\frac{\pi}{2} - \delta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}}}} & (20) \end{matrix}$

The second derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is:

$\begin{matrix} {\frac{d^{2}\left( {\Delta I_{\sum}^{2}} \right)}{{d\left( {❘{\Delta{\overset{.}{I}}_{n,y}}❘} \right)}^{2}} = {2 > 0}} & (21) \end{matrix}$

When Δİ*_(n,y) ends near the A₂ side, ΔI_(Σ.max) ² is:

$\begin{matrix} {{\Delta I_{\sum{.\max}}^{2}} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}} + {\Delta{\overset{.}{I}}_{n,y}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}^{2} + {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘}^{2} - {{2 \cdot {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}{\cos\left( {\frac{3\pi}{2} - \delta_{n}^{*}} \right)}}}}}} & (22) \end{matrix}$

At this moment, the derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is:

$\begin{matrix} {\frac{d\left( {\Delta I_{\sum}^{2}} \right)}{d\left( {❘{\Delta{\overset{.}{I}}_{n,y}}❘} \right)} = {{2{❘{\Delta{\overset{.}{I}}_{n,y}}❘}} - {2{{\cos\left( {\frac{3\pi}{2} - \delta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}}}} & (23) \end{matrix}$

When Δİ*_(n,y) ends near the A₂ side, the second derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is the same as formula (21):

Denoting Δİ_(n,y) ending at A₁ as Δİ_(n,y) ^(A1), and denoting Δİ_(n,y) ending at A₂ as Δİ_(n,y) ^(A2);

Respectively denoting derivatives in formula (20) and formula (23) as f_(n,y) ¹ and f_(n,y) ².

Preferably, step S7 of solving the optimization problem of the non-convex constraint space comprises the following specific steps:

Step S701: inputting data matrixes of İ_(1.tr)-İ_(n.tr), λ₁-λ_(n) and α₁-α_(n);

$\begin{matrix} \begin{bmatrix} I_{1.{tr}} & \beta_{1.{tr}} & \lambda_{1} & \alpha_{1} \\ I_{2.{tr}} & \beta_{2.{tr}} & \lambda_{2} & \alpha_{2} \\  \vdots & \vdots & \vdots & \vdots \\ I_{n.{tr}} & \beta_{n.{tr}} & \lambda_{n} & \alpha_{n} \end{bmatrix} & (28) \end{matrix}$

wherein I_(1.tr)-I_(n.tr) are amplitudes of İ_(1.tr)-İ_(n.tr) and β_(1.tr)-β_(1.tr) are phases of İ_(1.tr)-İ_(n.tr);

wherein the specific values of λ₁-λ_(n) and α₁-α_(n) can be determined by the type of a current transformer and the data synchronization accuracy of a data synchronization method, and İ_(1.tr)-İ_(n.tr) shall be fault current flowing through primary and secondary fusing intelligent switches when a fault outside the protected area is most severe;

Step S702: determining to discretize the constraint space of the optimization problem into (4+2m)n points according to λ₁-λ_(n) and α₁-α_(n);

Step S703: generating a discrete point matrix;

$\begin{matrix} \begin{bmatrix} {{I_{1.{tr}}\left( {1 - \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \alpha_{1}} \right)}} & {{I_{2.{tr}}\left( {1 - \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \alpha_{2}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 - \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \alpha_{n}} \right)}} \\ {{I_{1.{tr}}\left( {1 - \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \alpha_{1}} \right)}} & {{I_{2.{tr}}\left( {1 - \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \alpha_{2}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 - \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \alpha_{n}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \alpha_{1}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \alpha_{2}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \alpha_{n}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \alpha_{1}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \alpha_{2}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \alpha_{n}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \frac{\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \frac{\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \frac{\alpha_{n}}{2^{m}}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \frac{3\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \frac{3\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \frac{3\alpha_{n}}{2^{m}}} \right.}} \\  \vdots & \vdots & \ldots & \vdots \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \frac{\left( {{2k} + 1} \right)\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \frac{\left( {{2k} + 1} \right)\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \frac{\left( {{2k} + 1} \right)\alpha_{n}}{2^{m}}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \frac{\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \frac{\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \frac{\alpha_{n}}{2^{m}}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \frac{3\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \frac{3\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \frac{3\alpha_{n}}{2^{m}}} \right)}} \\  \vdots & \vdots & \ldots & \vdots \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \frac{\left( {{2k} + 1} \right)\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \frac{\left( {{2k} + 1} \right)\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \frac{\left( {{2k} + 1} \right)\alpha_{n}}{2^{m}}} \right)}} \end{bmatrix} & (29) \end{matrix}$

wherein k≥1 and 2k−1≤2^(m)−1

Step S704: for 1-n columns in the matrix in step S703, performing the following (4+2^(m))^(n) operations respectively for different points x_({dot over (ι)}1,1), x_({dot over (ι)}2,2), . . . , x_({dot over (ι)}n,n) of 4+2^(m) rows;

ΔI _(Σ) =|x _({dot over (ι)}1,1) +x _({dot over (ι)}2,2) + . . . +x _({dot over (ι)}n,n)|  (30)

wherein 1≤{dot over (ι)}1, {dot over (ι)}2, . . . ,{dot over (ι)}_(n)≤4+2^(m);

Step S705: searching for the maximum value in the (4+2^(m))^(n) results calculated in step S704, i.e., the approximate solution ΔI_(Σ.max)′ of a global optimal solution, and approximately replacing ΔI_(Σ.max) by ΔI_(Σ.max)′ as the protection setting value in formula (9).

Preferably, step S8 of analyzing the solution result of the optimization problem of the non-convex constraint space comprises:

Setting up a simulation model in PSCAD/EMTDC simulation software by using the system topology of a high-proportion distributed generator T-connected distribution network with a power electronic transformer as a transmission and distribution interface.

Preferably, step S8 of verifying the solution result of the optimization problem of the non-convex constraint space comprises:

Conducting contrast verification through results of sector boundary traversal and six-point, eight-point and twelve-point fast traversal methods.

Compared with the prior art, the present invention has the beneficial effects that:

The multi-terminal differential protection setting method for a distributed generator T-connected distribution network proposed by the preset invention can accurately identify internal faults and external faults of the protected area in the high-proportion distributed generator T-connected distribution network, with sufficient reliability and sensitivity, and can quickly and accurately calculate the overall maximum error of differential current of superimposed multi-terminal current phasor errors.

DESCRIPTION OF DRAWINGS

FIG. 1 is a schematic diagram of a high-proportion distributed generator T-connected distribution network of the present invention;

FIG. 2 is a schematic diagram of a configuration scheme of centralized multi-terminal current differential protection of the present invention;

FIG. 3 is a schematic diagram of the relationship between the measured value and the actual value of current of the present invention;

FIG. 4 is a schematic diagram of the measured value and the actual value of n-terminal current of the present invention;

FIG. 5 is a schematic diagram of the position relationship between Δİ*₁+Δİ*₂+ . . . +Δİ_(n-1) and Δİ*_(n) in a complex plane of the present invention;

FIG. 6 is a schematic diagram of numbering four boundaries of a sector shadow area corresponding to Δİ_(n) of the present invention;

FIG. 7 is a schematic diagram of reducing non-convex constraint space of an optimization problem from n sector shadow areas to 2n arcs in the present invention;

FIG. 8 is a schematic diagram of analysis of an internal arc of a sector passing through a connecting line A₁A₃ in the present invention;

FIG. 9 is a schematic diagram of reducing constraint space of an optimization problem from 2n arcs of a sector to n arcs and 2n points in the present invention;

FIG. 10 is a schematic diagram of drawing a tangent line of an arc through points A₂ and A₄ and the arc top by discretization of a boundary 4 using two characteristics in the present invention;

FIG. 11 is a schematic diagram of further reducing constraint space of an optimization problem from n arcs and 2n points to 6n points in the present invention;

FIG. 12 is a schematic diagram of drawing a tangent line of an arc through intersection points of OA₅ and OA₆ and the arc based on FIG. 11 in the present invention;

FIG. 13 is a schematic diagram of a high-proportion distributed generator T-connected distribution network with a power electronic transformer as a transmission and distribution interface of the present invention;

FIG. 14 is a flow chart of a multi-terminal differential protection setting method for a distributed generator T-connected distribution network of the present invention.

DETAILED DESCRIPTION

As shown in FIG. 14 , the present invention provides a multi-terminal differential protection setting method for a distributed generator T-connected distribution network, comprising the following steps:

S1: Connecting a distributed generator into a distribution network by means of T-connection, and making the distributed generator equivalent to a positive sequence current source controlled by grid interconnection voltage;

The distribution network of a high-proportion T-connected inverter interfaced distributed generator is shown in FIG. 1 . For the simplicity of analysis, only feeders to be studied are shown in FIG. 1 .

ESS is an energy storage station, PV is a distributed photovoltaic power station, LD is a T-connected load, and inverter interfaced distributed generators such as photovoltaic type and energy storage type are connected into the distribution network by means of T-connection.

An inverter interfaced distributed generator can be made equivalent to a positive sequence current source controlled by grid interconnection voltage under the low voltage ride-through control strategy and the negative sequence current suppression control strategy.

Further, in step S1, an equivalent model for making the distributed generator equivalent to a positive sequence current source controlled by grid interconnection voltage is:

$\begin{matrix} {{U_{{pcc}(1)}^{*} = {U_{{pcc}(1)}/U_{{{pcc}(1)}.N}}}{I_{{DGm}.N} = {\sqrt{2}P_{{ref}.N}/3U_{{{pcc}(1)}.N}}}{i_{q(1)} = \left\{ {{\begin{matrix} {0,} & {U_{{pcc}(1)}^{*} > 0.9} \\ {{{- 2}\left( {1 - U_{{pcc}(1)}^{*}} \right)I_{{DGm}.N}},} & {0.4 \leq U_{{pcc}(1)}^{*} \leq 0.9} \\ {{{- 1.2}I_{{DGm}.N}},} & {U_{{pcc}(1)}^{*} < 0.4} \end{matrix}i_{d(1)}} = {{{\min\left( {\frac{\sqrt{2}P_{ref}}{3U_{{pcc}(1)}},\sqrt{\left( {1.2I_{{DGm}.N}} \right)^{2} - i_{q(1)}^{2}}} \right)}{\overset{.}{I}}_{DG}} = \left\{ \begin{matrix} {{\frac{\sqrt{i_{d(1)}^{2} + i_{q(1)}^{2}}}{\sqrt{2}}{\angle\left( {{\arctan\frac{i_{q(1)}}{i_{d(1)}}} + \varphi} \right)}},} & {i_{d(1)} \neq 0} \\ {{\frac{\sqrt{i_{d(1)}^{2} + i_{q(1)}^{2}}}{\sqrt{2}}{\angle\left( {{{- 90}{^\circ}} + \varphi} \right)}},} & {i_{d(1)} = 0} \end{matrix} \right.}} \right.}} & (1) \end{matrix}$

wherein U_(pcc(1)) is an effective value of grid interconnection positive sequence phase voltage, U_(pcc(1).N) is an effective value of grid interconnection positive sequence phase voltage in the rated running state, I_(DGm.N) is an amplitude of rated current, P_(ref.N) is a rated active reference value, P_(ref) is an active reference value before failure, İ_(DG) is output current, {dot over (ι)}_(d(1)) and {dot over (ι)}_(q(1)) are a d-axis component and a q-axis component of the positive sequence component of the output current, and φ is the phase of {dot over (U)}_(pcc(1)).

It can be seen from formula (1) that the output current of the inverter interfaced distributed generator has the characteristics of limited amplitude and controlled phase. Under the action of a current limiting link, the amplitude of İ_(DG) will not exceed 1.2 times the rated value, either in normal operation or during system failure. Therefore, an inverter interfaced distributed generator connected into the distribution network by high-proportion T-connection has typical weak feed characteristics.

S2: Acquiring multi-terminal current according to the positive sequence current source controlled by grid interconnection voltage;

S3: Calculating the overall maximum error of superimposed multi-terminal current errors according to the multi-terminal current;

Aiming at the configuration scheme and action criterion of multi-terminal current differential protection of a high-proportion distributed generator T-connected distribution network, the present invention adopts the configuration scheme of centralized multi-terminal differential protection, as shown in FIG. 2 .

CB is a feeder outgoing circuit breaker, FSW1 to FSW3 are sectional switches on main feeders, and QSW1 to QSWn+1 are branch switches on T-connected lines, which are primary and secondary fusion intelligent switches. The protection system is composed of a master station, a plurality of primary and secondary fusion switches and a communication network. In order to ensure reliability, the network is generally constructed through optical fiber. After being synchronized with the data of the master station, the primary and secondary fusion switches complete data acquisition and upload electrical data to the master station, the master station completes the calculation of differential protection, decides whether each switch needs to trip, and then sends the trip signal to the primary and secondary fusion switches that need to trip, and the primary and secondary fusion switches act to clear the fault after receiving the trip signal.

Further, step S3 of calculating the overall maximum error of superimposed multi-terminal current errors according to the multi-terminal current comprises:

Expressing the power frequency current phasor involved in a group of n(n≥3)-terminal current differential protection by a set P:

P={İ ₁ ,İ ₂ , . . . ,İ _(n)}  (2)

The action criterion of the multi-terminal current differential protection is:

I _(op) ≥I _(set)  (3)

wherein I_(op) is differential current; and I_(set) is a protection setting value;

Expressing the relationship between the measured value and the actual value of each power frequency current phasor as follows:

İ _({dot over (ι)}.ms) =İ _({dot over (ι)}.tr) +Δİ _({dot over (ι)})  (4)

-   -   i=1,2, . . . ,n

wherein İ_({dot over (ι)}.ms) represents the measured value of the power frequency current phasor, İ_({dot over (ι)}.tr) represents the actual value of the power frequency current phasor, and Δİ_({dot over (ι)}) represents the error between the actual value and the measured value;

At this moment, the differential current is:

$\begin{matrix} \begin{matrix} {I_{op} = {❘{{\overset{.}{I}}_{1.{ms}} + {\overset{.}{I}}_{2.{ms}} + \ldots + {\overset{.}{I}}_{n.{ms}}}❘}} \\ {= {❘{{\overset{.}{I}}_{1.{tr}} + {\Delta{\overset{.}{I}}_{1}} + {\overset{.}{I}}_{2.{tr}} + {\Delta{\overset{.}{I}}_{2}} + \ldots + {\overset{.}{I}}_{n.{tr}} + {\Delta{\overset{.}{I}}_{n}}}❘}} \end{matrix} & (5) \end{matrix}$

According to the Kirchhoff s Current Law, meeting the formula below during normal operation and during failure outside a protected area:

İ _(1.tr) +İ _(2.tr) + . . . +İ _(n.tr)=0  (6)

According to formula (5) and formula (6), the differential current during normal operation and during failure outside a protected area is:

I _(op) =|Δİ ₁ +Δİ ₂ + . . . +Δİ _(n)|  (7)

Denoting the overall maximum error of superimposed multi-terminal current errors as ΔI_(Σ):

ΔI _(Σ) =|Δİ ₁ +Δİ ₂ + . . . +Δİ _(n)|  (8)

A protection setting value shall be selected to ensure that no maloperation of faults exist outside the protected area, and no refused operation of faults exist inside the protected area, which have a certain degree of distinction. Therefore, the protection setting value can be:

I _(set) =k _(rel) ΔI _(Σ.max)  (9)

wherein k_(rel) is a reliability coefficient, generally 1.1-1.3; and ΔI_(Σ.max) is the overall maximum error of superimposed multi-terminal current errors.

The overall maximum error ΔI_(Σ.max) of superimposed multi-terminal errors is related to multi-terminal current, and the specific acquisition method will be described in detail below.

The overall maximum error of differential current of superposed n-terminal current errors is analyzed as follows: assuming that under the influence of the measurement error and the data synchronization error, a power frequency current phasor has a maximum amplitude error of λ and a maximum phase error of α, the specific values are determined by the type of a current transformer and the data synchronization accuracy of a data synchronization method. In the case of errors, the relationship between İ_(ms) and İ_(tr) is shown in FIG. 3 .

wherein a solid line with an arrow represents an actual value, and a dotted line with an arrow represents a measured value. An arc {circumflex over (l)}₁ is symmetric about a straight line where İ_(tr) is located, with the head end O of İ_(tr) as the center of a circle, with (1−λ)I_(tr) as the radius and with 2α as the centering angle. An arc {circumflex over (l)}₂ is symmetric about a straight line where in is located, with O as the center of a circle, with (1+λ)I_(tr) as the radius and with 2α as the centering angle. A sector shadow area formed by {circumflex over (l)}₁ and {circumflex over (l)}₂ is a set of possible measured values.

The relationship between a measured value and an actual value can be described by the formula below:

$\begin{matrix} {{{1 - \lambda} \leq \frac{I_{ms}}{I_{tr}} \leq {1 + \lambda}}{{- \alpha} \leq {{{angle}\left( {\overset{.}{I}}_{ms} \right)} - {{angle}\left( {\overset{.}{I}}_{tr} \right)}} \leq \alpha}} & (10) \end{matrix}$

In the case of errors in n-terminal current phasors, as shown in FIG. 4 :

Assuming that when Δİ₁=Δİ*₁, Δİ₂=Δİ*₂, . . . , Δİ_(n-1)=, Δİ*_(n-1), Δİ_(n)=Δİ*_(n), the overall error ΔI_(Σ) reaches the maximum value ΔI_(Σ.max), that is:

ΔI _(Σ.max) =|Δİ* ₁ +Δİ* ₂ + . . . +Δİ* _(n-1) +Δİ* _(n)|  (11)

Obviously, solving ΔI*₁, ΔI*₂, . . . , ΔI*_(n-1), ΔI*_(n), and the overall maximum error ΔI_(Σ.max) is an optimization problem, and the objective function is:

ΔI _(Σ.max)=max(|Δİ ₁ +Δİ ₂ + . . . +Δİ _(n-1) +Δİ _(n)|)  (12)

The constraints of the optimization problem include:

$\begin{matrix} {{{1 - \lambda_{i}} \leq \frac{❘{{\overset{.}{I}}_{i.{tr}} + {\Delta{\overset{.}{I}}_{i}}}❘}{I_{i.{tr}}} \leq {1 + \lambda_{i}}}{{- \alpha_{i}} \leq {{{angle}\left( {{\overset{.}{I}}_{i.{tr}} + {\Delta{\overset{.}{I}}_{i}}} \right)} - {{angle}\left( {\overset{.}{I}}_{i.{tr}} \right)}} \leq \alpha_{i}}{{i = 1},2,\ldots,{n - 1},{n.}}} & (13) \end{matrix}$

It can be seen that the objective function of the optimization problem can be transformed into a quadratic function, but the inequality constraints of the optimization problem contain a large number of nonlinear functions, which is different from a general quadratic programming problem. Therefore, it is difficult to solve the optimization problem effectively by using conventional quadratic programming methods such as Lagrange method, Lemke method and interior point method.

Generally speaking, an iterative method shall be adopted to find the extreme values of a complex function under complex constraints. However, the inequality constraints of the optimization problem are corresponding to the sector shadow area in FIG. 4 , which is obviously non-convex space, so iterative algorithms applicable to convex constraint space, such as gradient descent method, steepest descent method and Newton method, are also difficult to play a role in the optimization problem.

For the optimization problem on the non-convex constraint space, a feasible method is to discretize the non-convex constraint space first, and then to solve the discretized constraint space. In order to minimize the existence range of the optimal solution of the non-convex constraint space, two characteristics of the non-convex space are analyzed below.

S4: Establishing non-convex constraint space according to the overall maximum error of superimposed multi-terminal current errors;

Further, step S4 of establishing non-convex constraint space according to the overall maximum error of superimposed multi-terminal current errors comprises: two characteristics of the non-convex constraint space;

First characteristic: the phasors of Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n) end on the boundaries of the respective sector shadow area.

With Δİ*_(n) as an example, assuming that position relationship between Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n) in a complex plane is show in FIG. 5 .

wherein θ*_(n) is the phase difference between phasors Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n), and θ*_(n) ∈[0, 2π).

According to the Law of Cosines, ΔI_(Σ.max) ² is:

$\begin{matrix} {{\Delta I_{\sum{.\max}}^{2}} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}}}❘}^{2} + {❘{\Delta{\overset{.}{I}}_{n}^{*}}❘}^{2} - {{2 \cdot {❘{\Delta{\overset{.}{I}}_{n}^{*}}❘} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}}}❘}}{\cos\left( {\pi - \theta_{n}^{*}} \right)}}}}} & (14) \end{matrix}$

Denoting Δİ_(n) a having the phase difference of θ*_(n) from Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) as Δİ_(n) ^(θ) ^(n) ^(*), and denoting Δİ_(n) ^(θ) ^(n) ^(*) ending on the boundaries of the shadow area as Δİ_(n) ^(θ) ^(n) ^(*,max).

Since Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ_(n) ^(θ) ^(n) ^(*) are independent of each other, the derivative of ΔI_(Σ) ² with respect to |Δİ_(n) ^(θ) ^(n) ^(*)| is:

$\begin{matrix} {\frac{d\left( {\Delta I_{\sum}^{2}} \right)}{d\left( {❘{\Delta{\overset{.}{I}}_{n}^{\theta_{n}^{*}}}❘} \right)} = {{2{❘{\Delta{\overset{.}{I}}_{n}^{\theta_{n}^{*}}}❘}} - {2 \cdot {\cos\left( {\pi - \theta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}}}❘}}}} & (15) \end{matrix}$

The second derivative of ΔI_(Σ) ² with respect to |Δİ_(n) ^(θ) ^(n) ^(*)| is:

$\begin{matrix} {\frac{d^{2}\left( {\Delta I_{\sum}^{2}} \right)}{{d\left( {❘{\Delta{\overset{.}{I}}_{n}^{\theta_{n}^{*}}}❘} \right)}^{2}} = {2 > 0}} & (16) \end{matrix}$

According to formula (14) to formula (16), since the quadratic term coefficients of the second derivative and the first derivative are both larger than 0, as |Δİ_(n) ^(θ) ^(n) ^(*)| changes from 0 to |Δİ_(n) ^(θ) ^(n) ^(*,max)|, ΔI_(Σ) has at most one extreme point which is a minimum point. Therefore, ΔI_(Σ) is maximum only when Δİ_(n) ^(θ) ^(n) ^(*)=0 or Δİ_(n) ^(θ) ^(n) ^(*)=Δİ_(n) ^(θ) ^(n) ^(*,max).

Since Δİ_(n) which has the same direction as Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and whose amplitude is not zero can always be found, at least one Δİ_(n) ^(ξ) inevitably exists in the sector shadow area so that:

|Δİ* ₁ +Δİ* ₂ + . . . +Δİ* _(n-1) +Δİ _(n) ^(ξ) |>|Δİ* ₁ +Δİ ₂ + . . . +Δİ* _(n-1)+0|  (17)

Therefore, Δİ*_(n) n at the overall maximum error will not be 0, and thus can only be Δİ*_(n)=Δİ_(n) ^(θ) ^(n) ^(*,max).

In conclusion, when the overall error is maximum, the phasor of Δİ*_(n) must end on the boundaries of the sector shadow area. Obviously, when Δİ*₁-Δİ*_(n-1) are analyzed respectively, the same conclusion as Δİ*_(n) also can be made. Therefore, when the overall error is maximum, the phasors of Δİ*₁-Δİ*_(n) end on the boundaries of the sector shadow area.

Second characteristic: the phasors of Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n) must be on both ends of the straight boundaries if ending on the straight boundaries of the sector shadow area.

With Δİ*_(n) as an example, the second characteristic will be proved below.

Numbering four boundaries of the sector shadow area corresponding to Δİ_(n), as shown in FIG. 6 .

In the figure, the boundaries 1 and 2 are straight boundaries, and the boundaries 3 and 4 are arc boundaries, wherein Δİ_(n) ending on the boundary 1 can be decomposed into:

Δİ _(n) =Δİ* _(n,x) +Δİ _(n,y)  (18)

wherein Δİ*_(n,x) is the projection of Δİ_(n) on the vertical line of the boundary 1, and Δİ_(n,y) is the projection of Δİ_(n) on the boundary 1.

When Δİ*_(n,y) ends near the A₁ side, ΔI_(Σ.max) ² is:

$\begin{matrix} {{\Delta I_{\sum{.\max}}^{2}} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}} + {\Delta{\overset{.}{I}}_{n,y}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}^{2} + {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘}^{2} - {{2 \cdot {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}{\cos\left( {\frac{\pi}{2} - \delta_{n}^{*}} \right)}}}}}} & (19) \end{matrix}$

wherein δ*_(n) is the phase difference between Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n,x), and δ*_(n)∈[0, 2π).

Since Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n,x) Δİ_(n,y) are independent of each other, the derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is:

$\begin{matrix} {\frac{d\left( {\Delta I_{\sum}^{2}} \right)}{d\left( {❘{\Delta{\overset{.}{I}}_{n,y}}❘} \right)} = {{2{❘{\Delta{\overset{.}{I}}_{n,y}}❘}} - {2{{\cos\left( {\frac{\pi}{2} - \delta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}}}} & (20) \end{matrix}$

The second derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is:

$\begin{matrix} {\frac{d^{2}\left( {\Delta I_{\sum}^{2}} \right)}{{d\left( {❘{\Delta{\overset{.}{I}}_{n,y}}❘} \right)}^{2}} = {2 > 0}} & (21) \end{matrix}$

When Δİ*_(n,y) ends near the A₂ side, ΔI_(Σ.max) ² is:

$\begin{matrix} {{\Delta I_{\sum{.\max}}^{2}} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}} + {\Delta{\overset{.}{I}}_{n,y}^{*}}}❘}^{2} = {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}^{2} + {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘}^{2} - {{2 \cdot {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}{\cos\left( {\frac{3\pi}{2} - \delta_{n}^{*}} \right)}}}}}} & (22) \end{matrix}$

At this moment, the derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is:

$\begin{matrix} {\frac{d\left( {\Delta I_{\sum}^{2}} \right)}{d\left( {❘{\Delta{\overset{.}{I}}_{n,y}}❘} \right)} = {{2{❘{\Delta{\overset{.}{I}}_{n,y}}❘}} - {2{{\cos\left( {\frac{3\pi}{2} - \delta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}}}} & (23) \end{matrix}$

When Δİ*_(n,y) ends near the A₂ side, the second derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is the same as formula (21).

Denoting Δİ_(n,y) ending at A₁ as Δİ_(n,y) ^(A1), and denoting Δİ_(n,y) ending at A₂ as Δİ_(n,y) ^(A2).

Respectively denoting derivatives in formula (20) and formula (23) as f_(n,y) ¹ and f_(n,y) ².

According to formula (19) to formula (23), since the quadratic term coefficients of the second derivative and the first derivative are both larger than 0, as |Δİ_(n,y)|, changes from 0 to |Δİ_(n,y) ^(A1)|, ΔI_(Σ) has at most one extreme point which is a minimum point. As |Δİ_(n,y)| changes from 0 to |Δİ_(n,y) ^(A2)|, ΔI_(Σ) also has at most one extreme point which is also a minimum point. Therefore, ΔI_(Σ) is maximum only when Δİ_(n,y)=0, Δİ_(n,y)=Δİ_(n,y) ^(A1) or Δİ_(n,y)=Δİ_(n,y) ^(A2).

However, the formula below must be true regardless of the value of δ*_(n):

$\begin{matrix} {{{\cos\left( {\frac{3\pi}{2} - \delta_{n}^{*}} \right)} \cdot {\cos\left( {\frac{\pi}{2} - \delta_{n}^{*}} \right)}} \leq 0} & (24) \end{matrix}$

Therefore, at least one of f_(n,y) ¹ and f_(n,y) ² is larger than or equal to 0 regardless of the value of |Δİ_(n,y)|, that is, the formula below is true.

$\begin{matrix} {f_{n,y}^{1} = {{{2{❘{\Delta{\overset{.}{I}}_{n,y}}❘}} - {2{{\cos\left( {\frac{\pi}{2} - \delta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}}} \geq {0{or}}}} & (25) \end{matrix}$ $f_{n,y}^{2} = {{{2{❘{\Delta{\overset{.}{I}}_{n,y}}❘}} - {2{{\cos\left( {\frac{3\pi}{2} - \delta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}}} \geq 0}$

Therefore, when Δİ_(n,y)=Δİ_(n,y) ^(A1) or Δİ_(n,y)=Δİ_(n,y) ^(A2), the value of ΔI_(Σ) is larger than the value of ΔI_(Σ) when Δİ_(n,y)=0. That is, when the overall error is maximum, the phasor of ΔI*_(n) must be on both ends of the straight boundaries if ending on the straight boundaries of the sector shadow area.

Obviously, when Δİ*₁-Δİ*_(n-1) are analyzed respectively, the same conclusion as Δİ*_(n) also can be made. Therefore, the phasors of Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n) must be on both ends of the straight boundaries if ending on the straight boundaries of the sector shadow area.

S5: Reducing and discretizing the non-convex constraint space;

S6: Obtaining the optimization problem of the non-convex constraint space based on the reducing and discretizing results of the non-convex constraint space;

According to the two characteristics, the non-convex constraint space of the optimization problem can be reduced from n sector shadow areas to 2n arcs, as shown in FIG. 7 .

At this moment, the range of the optimal solution is still a little big, and further simplification is needed for the two arcs of the sector.

With Δİ*_(n) as an example, the internal arc of the sector, i.e., the boundary 3, can be analyzed through the connecting line A₁A₃, as shown in FIG. 8 .

In the figure, a is any point on the boundary 3, and Δİ_(n) ending at a is extended to _(i)n_(t)ersect with the line A₁A₃ at b.

According to the first characteristic, since any point on the line A₁A₃ is outside the sector area, the overall error ΔI_(Σ) when Δİ_(n) ends at b is larger than that when Δİ_(n) ends at a. According to the second characteristic, it can be concluded that the overall error ΔI_(Σ) when Δİ_(n) ends at A₁ or A₃ is larger than that when n ends at b.

Therefore, the constraint space of the optimization problem can be reduced from 2n arcs of the sector to n arcs and 2n points, as shown in FIG. 9 .

For the external arc of the sector, i.e., the boundary 4, since any point on the line A₂A₄ is inside the sector area, the boundary 4 cannot be simply reduced to two endpoints of the arc according to the two characteristics. To discretize the boundary 4 using the two characteristics, a tangent line of the arc can be drawn respectively through points A₂ and A₄ and the arc top, as shown in FIG. 10 .

In the figure, a tangent line through the point A₂ and a tangent line through the arc top A₇ intersect at A₅, and a tangent line through the point A₄ and a tangent line through the arc top A₇ intersect at A₆. It is not difficult to show that both OA₅ and OA₆ are angular bisectors of an angle α.

Obviously, any point on the line A₂A₅, the line A₅A₆ and the line A₄A₆ is outside the sector area, so the boundary 4 can be discretized into four points A₂, A₅, A₆ and A₄ according to the two characteristics.

Therefore, the constraint space of the optimization problem can be further reduced from n arcs and 2n points to 6n points, as shown in FIG. 11 .

However, as a result, the range of the constraint space is expanded, the solution obtained is only the approximate solution ΔI_(Σ.max)′ of the global optimal solution, but not a true global optimal solution ΔI_(Σ.max), and the value of ΔI_(Σ.max)′ may be larger than ΔI_(Σ.max). The difference between ΔI_(Σ.max)′ and ΔI_(Σ.max) is closely related to the distance from points on the line A₂A₅, the line A₅A₆ and the line A₄A₆ to the arc. In FIG. 10 , the farthest points from the arc are A₅ and A₆, and with I_(tr) as the reference value, the per unit distance from the points A₅ and A₆ to the arc is:

$\begin{matrix} {d = {\left( {1 + \lambda} \right)\left( {\frac{1}{\cos\left( {\alpha/2} \right)} - 1} \right)}} & (26) \end{matrix}$

It can be seen that the smaller λ and α are, the smaller d is. For example, when λ =0.1 and α=5°, d=0.001047, the value is small, then ΔI_(Σ.max)′ and ΔI_(Σ.max) are almost the same, the difference therebetween is within an acceptable range. However, as λ and α become larger, the difference between ΔI_(Σ.max)′ and ΔI_(Σ.max) will become larger. To solve the problem, a tangent line of the arc can be drawn respectively through intersection points of OAs and OA₆ and the arc based on FIG. 11 , as shown in FIG. 12 .

In the figure, A₈, A₉, A₁₀ and A₁₁ are intersection points of new tangent lines.

At this moment, the farthest points from the arc are A₈ to A₁₁, and with I_(tr) as the reference value, the per unit distance is:

$\begin{matrix} {d = {\left( {1 + \lambda} \right)\left( {\frac{1}{\cos\left( {\alpha/4} \right)} - 1} \right)}} & (27) \end{matrix}$

Compared with formula (26), the distance d decreases, and the distance between ΔI_(Σ.max)′ and ΔI_(Σ.max) will also decrease. As λ and α become larger, a similar method can be adopted to continuously draw tangents to select intersection points of the tangent lines, the arc is replaced by multiple tangent lines, i.e., the curved line is replaced by straight lines, and then the boundary 4 can be discretized into intersection points of tangent lines by the two characteristics. As a result, the constraint space of the optimization problem can be discretized into 6n, 8n, 12n, 20n, . . . , (4+2m)n points, wherein m≥1. The more the points are, the closer the approximate solution of the global optimal solution is to the true global optimal solution. However, no matter how many points are discretized, the approximate solution of the global optimal solution will always be slightly larger than the true global optimal solution. Therefore, even if the true global optimal solution cannot be obtained, the approximate solution of the global optimal solution can meet the need of protection setting.

S7: Solving the optimization problem of the non-convex constraint space;

Further, step S7 of solving the optimization problem of the non-convex constraint space comprises the following specific steps:

Step S701: inputting data matrixes of İ_(1.tr)-İ_(n.tr), λ₁-λ_(n) and α₁-α_(n);

$\begin{matrix} \begin{bmatrix} I_{1.{tr}} & \beta_{1,{tr}} & \lambda_{1} & \alpha_{1} \\ I_{2.{tr}} & \beta_{2,{tr}} & \lambda_{2} & \alpha_{2} \\  \vdots & \vdots & \vdots & \vdots \\ I_{n,{tr}} & \beta_{n,{tr}} & \lambda_{n} & \alpha_{n} \end{bmatrix} & (28) \end{matrix}$

wherein I_(1.tr)-I_(n.tr) are amplitudes of İ_(1.tr)-İ_(n.tr), and β_(1.tr)-β_(1.tr) are phases of İ_(1.tr)-İ_(n.tr);

wherein the specific values of λ₁-λ_(n) and α₁-α_(n) can be determined by the type of a current transformer and the data synchronization accuracy of a data synchronization method, and İ_(1.tr)-İ_(n.tr) shall be fault current flowing through primary and secondary fusing intelligent switches when a fault outside the protected area is most severe;

Step S702: determining to discretize the constraint space of the optimization problem into (4+2m)n points according to λ₁-λ_(n) and α₁-α_(n);

Step S703: generating a discrete point matrix;

$\begin{matrix} \begin{bmatrix} {{I_{1.{tr}}\left( {1 - \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \alpha_{1}} \right)}} & {I_{2.{tr}}\left( {1 - \lambda_{2}} \right)\angle\left( {\beta_{2.{tr}} - \alpha_{2}} \right)} & \ldots & {I_{n.{tr}}\left( {1 - \lambda_{n}} \right)\angle\left( {\beta_{n.{tr}} - \alpha_{n}} \right)} \\ {I_{1.{tr}}\left( {1 - \lambda_{1}} \right)\angle\left( {\beta_{1.{tr}} + \alpha_{1}} \right)} & {I_{2.{tr}}\left( {1 - \lambda_{2}} \right)\angle\left( {\beta_{2.{tr}} + \alpha_{2}} \right)} & \ldots & {I_{n.{tr}}\left( {1 - \lambda_{n}} \right)\angle\left( {\beta_{n.{tr}} + \alpha_{n}} \right)} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \alpha_{1}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \alpha_{2}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \alpha_{n}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \alpha_{1}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \alpha_{2}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \alpha_{n}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \frac{\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \frac{\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \frac{\alpha_{n}}{2^{m}}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \frac{3\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \frac{3\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \frac{3\alpha_{n}}{2^{m}}} \right)}} \\  \vdots & \vdots & \ldots & \vdots \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \frac{\left( {{2k} + 1} \right)\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \frac{\left( {{2k} + 1} \right)\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \frac{\left( {{2k} + 1} \right)\alpha_{n}}{2^{m}}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \frac{\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \frac{\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \frac{\alpha_{n}}{2^{m}}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \frac{3\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \frac{3\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \frac{3\alpha_{n}}{2^{m}}} \right)}} \\  \vdots & \vdots & \ldots & \vdots \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \frac{\left( {{2k} + 1} \right)\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \frac{\left( {{2k} + 1} \right)\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \frac{\left( {{2k} + 1} \right)\alpha_{n}}{2^{m}}} \right)}} \end{bmatrix} & (29) \end{matrix}$

wherein k≥1 and 2k−1≤2^(m)−1

Step S704: for 1-n columns in the matrix in step S703, performing the following (4+2^(m))^(n) operations respectively for different points x_({dot over (ι)}1,1), x_({dot over (ι)}2,2), . . . , x_({dot over (ι)}n,n) of 4+2^(m) rows;

ΔI _(Σ) =|x _({dot over (ι)}1,1) +x _({dot over (ι)}2,2) + . . . +x _({dot over (ι)}n,n)|  (30)

wherein 1≤{dot over (ι)}1, {dot over (ι)}2, . . . , {dot over (ι)}_(n)≤4+2^(m);

Step S705: searching for the maximum value in the (4+2^(m))^(n) results calculated in step S704, i.e., the approximate solution ΔI_(Σ.max)′ of a global optimal solution, and approximately replacing ΔI_(Σ.max) by ΔI_(Σ.max)′ as the protection setting value in formula (9).

S8: Analyzing and verifying the solution result of the optimization problem of the non-convex constraint space.

The power electronic transformer is considered as ideal equipment to replace a conventional transformer due to significant advantages in terms of power allocation and voltage regulation, which can realize flexible regulation of power resources from the distribution network to transmission networks as well as active absorption and controllable transmission of high-proportion new energy in the distribution network. Due to the limited discharge capacity, the power electronic transformer will also have the function of limiting fault current in order to ensure no burning during failure. For the high-proportion distributed generator T-connected distribution network with a power electronic transformer as a transmission and distribution interface, all the power supplies in the system will show typical weak feed characteristics after failure, forming a fully weak-fed distribution network. The fully weak-fed distribution network is a new type of high-proportion distributed generator T-connected distribution network, which has more stringent requirements for protection setting than the distribution network with a conventional transformer as a transmission and distribution interface. Therefore, the high-distributed generator T-connected distribution network with a power electronic transformer as a transmission and distribution interface is taken as an example for analysis.

A simulation model is set up in PSCAD/EMTDC simulation software by using the system topology of the high-proportion distributed generator T-connected distribution network with a power electronic transformer as a transmission and distribution interface according to FIG. 13 .

The distribution network has a voltage class of 10 kV and a rated frequency of 50 Hz, and the system adopts a neutral ungrounding mode. The power electronic transformer has a capacity of 10 MW, adopts the strategy of constant voltage constant frequency control/voltage dependent current limit control, and limits the fault current amplitude to less than 2 times the rated current when a short circuit fault occurs. PV1-PV3 and ESS1-ESS2 all have a rated capacity of 1 MW, output power of 1 p.u. per unit power factor before failure, and adopt the PQ control strategy, the negative sequence current suppression control strategy and the low voltage ride-through control strategy considering fault current limiting, with current limiting of 1.2 p.u. The impedance of the line is Z₁=0.1+j0.3 Ω/km and Z₀=0.3+j0.9 Ω/km, the impedance of T-connected lines is 1 km, and the impedance interval between adjacent T-connected lines is 2 km; and the impedance of LD1-LD3 is 2 MW. With a group of multi-terminal differential protection composed of FSW1, FSW2 and QSW2-QSW4 as an example, f1-f5 have faults outside the protected area, and f6 has faults inside the protected area.

The effectiveness of the proposed method for solving the overall maximum error of superposed multi-terminal errors is verified. For three-phase grounded short circuit at f1-f5, the current simulation values of intelligent switches are shown in Table 1, with the unit of kA (which are A-phase current, with the A-phase current phase of FSW1 as reference).

TABLE 1 Current Simulation Value of Intelligent Switches in the Case of Three-Phase Short Circuit at Different Positions FSW1 QSW2 QSW3 QSW4 FSW2 f1 0.2714∠0.0º 0.0689∠192.9º 0.0693∠184.3º 0.0684∠172.6º 0.0684∠170.1º f2 1.1028∠0.0º 0.9001∠178.3º 0.0693∠201.9º 0.0695∠183.1º 0.0695∠178.0º f3 1.0405∠0.0º 0.0709∠−41.1º 1.0053∠171.8º 0.0694∠227.6º 0.0694∠221.2º f4 0.9803∠0.0º 0.0705∠−49.1º 0.0693∠−31.0º 1.0675∠171.7º 0.0694∠242.2º f5 0.9288∠0.0º 0.0809∠−58.2º 0.0693∠−29.1º 0.0710∠−47.3º 1.0910∠171.9º

Assuming that all the intelligent switches have a maximum amplitude error of λ=10% and a maximum angle error of α=5°, when three-phase grounded short circuit occurs at f1-f5, the overall maximum error of differential current of superimposed multi-terminal current errors calculated by sector area traversal, sector boundary traversal and the six-point, eight-point and twelve-point fast traversal methods proposed by the present invention is shown in Table 2, wherein the sector area traversal means, in the whole sector area, selecting one point at intervals of 2% of amplitude and 1° of angle, 121 points in total; the sector boundary traversal means, on the sector boundary, selecting one point at intervals of 2% of amplitude on the straight boundary and of 1° on the arc boundary, 40 points in total; and the six-point, eight-point and twelve-point fast traversal methods are calculation methods described in the embodiment. With a computer with an Intel® Core™ i7-6700HQ CPU @2.60 GHz 2.59 GHz processor, an algorithm program is written by matlab2021b commercial software package for calculation.

TABLE 2 Overall Maximum Error and Calculation Time of Superposed Multi- terminal Errors under Different Calculation Methods in the Case of Three-Phase Short Circuit at Different Positions Full Boundary Twelve- Eight- Six- Traversal Traversal Point Point Point f1 Maximum error 0.0741 0.0741 0.0741 0.0741 0.0741 (kA) Calculation 3925.053 26.940 0.141 0.063 0.058 time (s) f2 Maximum error 0.2997 0.2997 0.2997 0.2997 0.2997 (kA) Calculation 3852.813 27.138 0.157 0.076 0.051 time (s) f3 Maximum error 0.3507 0.3507 0.3507 0.3507 0.3507 (kA) Calculation 4012.235 26.993 0.147 0.063 0.052 time (s) f4 Maximum error 0.3065 0.3065 0.3065 0.3065 0.3065 (kA) Calculation 3965.398 26.326 0.136 0.061 0.050 time (s) f5 Maximum error 0.3037 0.3037 0.3037 0.3037 0.3037 (kA) Calculation 3896.897 26.856 0.142 0.067 0.056 time (s)

By comparing the results of sector area traversal with those of sector boundary traversal, it can be seen that the first characteristic introduced in Part 3 of the present invention is correct. By comparing the results of sector boundary traversal with those of the six-point, eight-point and twelve-point fast traversal methods, the correctness of the second characteristic introduced in Part 3 and the effectiveness of the proposed fast traversal calculation methods based on high discretization of non-convex space are verified. By comparing the calculation time of different calculation methods, it can be seen that the calculation time of the proposed fast traversal calculation methods is shorter, which is 4-5 orders of magnitude less than that of the sector area traversal method and the sector boundary traversal method. In order to give consideration to the calculation speed and the accuracy of the calculation results, the twelve-point fast traversal calculation method can be used for protection setting.

In order to ensure the reliability of the proposed protection scheme without misoperation in case of faults outside the protected area, the maximum value of the overall maximum error of superposed multi-terminal errors in the case pf serious three-phase grounded short circuit faults at different positions shall be selected, i.e., 0.3065 kA, multiplied by a reliability coefficient k_(rel)=1.1, as the protection setting value I_(set)=1.1×0.3065=0.3372 kA of the proposed protection scheme.

In order to analyze the reliability of the proposed protection scheme without misoperation in case of faults outside the area, a three-phase short circuit fault is set at f1-f5 outside the protected area in the simulation model. On the basis of the fault current simulation value of each intelligent switch, an amplitude error of 10% and an angle error of 5° are superimposed, so as to obtain the fault current measured value of each intelligent switch, and on the basis of the measured value, the difference current of each phase is calculated and compared with the protection setting value to determine the action, as shown in Table 3.

TABLE 3 Differential Current and Protection Action of Three- Phase Short Circuit Fault outside Protected Area f1 f2 f3 f4 f5 Differential Phase A 0.0546 0.2211 0.2256 0.2370 0.2241 current Phase B 0.0476 0.1927 0.1966 0.1968 0.1953 (kA) Phase C 0.0511 0.2070 0.2112 0.2163 0.2117 Protection No No No No No action action action action action action

It can be seen that when each intelligent switch has an amplitude error of 10% and an angle error of 5°, the protection will not misoperate when a serious three-phase short circuit fault occurs outside the protected area, so the proposed protection scheme has good reliability.

In order to analyze the sensitivity of the proposed protection scheme when a fault occurs inside the protected area, a phase-to-phase fault of the two phases A and B is set at f6 inside the protected area of the simulation model, and the transition resistance Z_(f) is respectively 1Ω, 5Ω, 10Ω and 20Ω. On the basis of the fault current simulation value of each intelligent switch, an amplitude error of 10% and an angle error of 5° are superimposed to obtain the fault current measured value of each intelligent switch, and on the basis of the measured value, the difference current of each phase is calculated and compared with the protection setting value to determine the action, as shown in Table 4.

TABLE 4 Differential Current and Protection Action of Phase-To-Phase Fault of Two Phases inside Protected Area Z_(f) = 1Ω Z_(f) = 5Ω Z_(f) = 10Ω Z_(f) = 20Ω Differential current Phase A 0.8554 0.7276 0.5288 0.3860 (kA) Phase B 0.8493 0.7114 0.5202 0.3749 Phase C 0.0704 0.0649 0.0603 0.0523 Protection action Action Action Action Action

It can be seen that when each intelligent switch has an amplitude error of 10% and an angle error of 5°, the protection can act reliably when a phase-to-phase short circuit fault with larger transition resistance occurs in the protected area, and the proposed protection scheme has good sensitivity.

In addition, the conditions with an amplitude error less than 10% and an angular error less than 5° are randomly superposed in the simulation values of different fault positions and different fault types as the measured value to conduct a large amount of example analysis. The results show that the multi-terminal differential protection setting scheme for a distribution network proposed by the present invention has sufficient sensitivity and reliability.

The present invention proposes a fast traversal method based on high discretization of non-convex space, which can quickly and accurately calculate the overall maximum error of differential current of superimposed multi-terminal current phasor errors, settling the problem of solving the setting value of multi-terminal current differential protection.

The multi-terminal differential protection setting method for a distributed generator T-connected distribution network proposed by the preset invention can accurately identify internal faults and external faults of the protected area in the high-proportion distributed generator T-connected distribution network, with sufficient reliability and sensitivity, and can quickly and accurately calculate the overall maximum error of differential current of superimposed multi-terminal current phasor errors. 

What is claimed is:
 1. A multi-terminal differential protection setting method for a distributed generator T-connected distribution network, comprising the following steps: S1: connecting a distributed generator into a distribution network by means of T-connection, and making the distributed generator equivalent to a positive sequence current source controlled by grid interconnection voltage; S2: acquiring multi-terminal current according to the positive sequence current source controlled by grid interconnection voltage; S3: calculating the overall maximum error of superimposed multi-terminal current errors according to the multi-terminal current; S4: establishing non-convex constraint space according to the overall maximum error of superimposed multi-terminal current errors; S5: reducing and discretizing the non-convex constraint space; S6: obtaining the optimization problem of the non-convex constraint space based on the reducing and discretizing results of the non-convex constraint space; S7: solving the optimization problem of the non-convex constraint space; S8: analyzing and verifying the solution result of the optimization problem of the non-convex constraint space.
 2. The multi-terminal differential protection setting method for a distributed generator T-connected distribution network according to claim 1, wherein in step S1, an equivalent model for making the distributed generator equivalent to a positive sequence current source controlled by grid interconnection voltage is: $\begin{matrix} {U_{{pcc}(1)}^{*} = {U_{{pcc}(1)}/U_{{{pcc}(1)}.N}}} & (1) \end{matrix}$ $I_{{DGm}.N} = {\sqrt{2}P_{{ref}.N}/3U_{{{pcc}(1)}.N}}$ $i_{q(1)} = \left\{ \begin{matrix} {0,} & {U_{{pcc}(1)}^{*} > 0.9} \\ {{{- 2}\left( {1 - U_{{pcc}(1)}^{*}} \right)I_{{DGm}.N}},} & {0.4 \leq U_{{pcc}(1)}^{*} \leq 0.9} \\ {{{- 1.2}I_{{DGm}.N}},} & {U_{{pcc}(1)}^{*} < 0.4} \end{matrix} \right.$ $i_{d(1)} = {\min\left( {\frac{\sqrt{2}P_{ref}}{3U_{{pcc}(1)}},\sqrt{\left. {\left( {1.2I_{{DGm}.N}} \right)^{2} - i_{q(1)}^{2}} \right)}} \right.}$ ${\overset{.}{I}}_{DG} = \left\{ \begin{matrix} {{\frac{\sqrt{i_{d(1)}^{2} + i_{q(1)}^{2}}}{\sqrt{2}}{\angle\left( {{\arctan\frac{i_{q(1)}}{i_{d(1)}}} + \varphi} \right)}},} & {i_{d(1)} \neq 0} \\ {{\frac{\sqrt{i_{d(1)}^{2} + i_{q(1)}^{2}}}{\sqrt{2}}{\angle\left( {{{- 90}{^\circ}} + \varphi} \right)}},} & {i_{d(1)} = 0} \end{matrix} \right.$ wherein U_(pcc(1)) is an effective value of grid interconnection positive sequence phase voltage, U_(pcc(1).N) is an effective value of grid interconnection positive sequence phase voltage in the rated running state, I_(DGm.N) is an amplitude of rated current, P_(ref.N) is a rated active reference value, P_(ref) is an active reference value before failure, İ_(DG) is output current, {dot over (ι)}_(d(1)) and {dot over (ι)}_(q(1)) are a d-axis component and a q-axis component of the positive sequence component of the output current, and φ is the phase of {dot over (U)}_(pcc(1)).
 3. The multi-terminal differential protection setting method for a distributed generator T-connected distribution network according to claim 1, wherein step S3 of calculating the overall maximum error of superimposed multi-terminal current errors according to the multi-terminal current comprises: expressing the power frequency current phasor involved in a group of n(n≥3)-terminal current differential protection by a set P: P={İ ₁ ,İ ₂ , . . . ,İ _(n)}  (2) the action criterion of the multi-terminal current differential protection is: I _(op) ≥I _(set)  (3) wherein I_(op) is differential current; and I_(set) is a protection setting value; expressing the relationship between the measured value and the actual value of each power frequency current phasor as follows: İ _({dot over (ι)}.ms) =İ _({dot over (ι)}.tr) +Δİ _({dot over (ι)})  (4) i=1,2, . . . , n wherein İ_({dot over (ι)}.ms) represents the measured value of the power frequency current phasor, İ_({dot over (ι)}.tr) represents the actual value of the power frequency current phasor, and Δİ_({dot over (ι)}) represents the error between the actual value and the measured value; at this moment, the differential current is: $\begin{matrix} \begin{matrix} {I_{op} = {❘{{\overset{.}{I}}_{1.{ms}} + {\overset{.}{I}}_{2.{ms}} + \ldots + {\overset{.}{I}}_{n.{ms}}}❘}} \\ {= {❘{{\overset{.}{I}}_{1.{tr}} + {\Delta{\overset{.}{I}}_{1}} + {\overset{.}{I}}_{2.{tr}} + {\Delta{\overset{.}{I}}_{2}} + \ldots + {\overset{.}{I}}_{n.{tr}} + {\Delta{\overset{.}{I}}_{n}}}❘}} \end{matrix} & (5) \end{matrix}$ according to the Kirchhoffs Current Law, meeting the formula below during normal operation and during failure outside a protected area: İ _(1.tr) +İ _(2.tr) + . . . +İ _(n.tr)=0  (6) according to formula (5) and formula (6), the differential current during normal operation and during failure outside a protected area is: I _(op) =|Δİ ₁ +Δİ ₂ + . . . +Δİ _(n)|  (7) denoting the overall maximum error of superimposed multi-terminal current errors as ΔI_(Σ): ΔI _(Σ) =|Δİ ₁ +Δİ ₂ + . . . +Δİ _(n)|  (8)
 4. The multi-terminal differential protection setting method for a distributed generator T-connected distribution network according to claim 1, wherein step S4 of establishing non-convex constraint space according to the overall maximum error of superimposed multi-terminal current errors comprises: according to Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n), establishing a sector shadow area, wherein the phasors of Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n) end on the boundaries of the respective sector shadow area; if the position relationship between Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n) in a complex plane is that Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n) are at an angle of θ*_(n); wherein θ*_(n) is the phase difference between phasors Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n), and θ*_(n)∈[0, 2π); according to the Law of Cosines, ΔI_(Σ.max) ² is: $\begin{matrix} \begin{matrix} {{\Delta I_{\Sigma.\max}^{2}} = {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n}^{*}}}❘}^{2}} \\ {= {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}}}❘}^{2} + {❘{\Delta{\overset{.}{I}}_{n}^{*}}❘}^{2} - {2 \cdot}}} \\ {{{❘{\Delta{\overset{.}{I}}_{n}^{*}}❘} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}}}❘}}{\cos\left( {\pi - \theta_{n}^{*}} \right)}} \end{matrix} & (14) \end{matrix}$ denoting Δİ_(n) having the phase difference of θ*_(n) from Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) as Δİ_(n) ^(θ) ^(n) ^(*), and denoting Δİ_(n) ^(θ) ^(n) ^(*) ending on the boundaries of the shadow area as Δİ_(n) ^(θ) ^(n) ^(*, max); since Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ_(n) ^(θ) ^(n) ^(*) are independent of each other, the derivative of ΔI_(Σ) ² with respect to |Δİ_(n) ^(θ) ^(n) ^(*)| is: $\begin{matrix} {\frac{d\left( {\Delta I_{\Sigma}^{2}} \right)}{d\left( {❘{\Delta{\overset{.}{I}}_{n}^{\theta_{n}^{*}}}❘} \right)} = {{2{❘{\Delta{\overset{.}{I}}_{n}^{\theta_{n}^{*}}}❘}} - {2 \cdot {\cos\left( {\pi - \theta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}}}❘}}}} & (15) \end{matrix}$ the second derivative of ΔI_(Σ) ² with respect to |Δİ_(n) ^(θ) ^(n) ^(*)| is: $\begin{matrix} {\frac{d^{2}\left( {\Delta I_{\Sigma}^{2}} \right)}{{d\left( {❘{\Delta{\overset{.}{I}}_{n}^{\theta_{n}^{*}}}❘} \right)}^{2}} = {2 > 0}} & (16) \end{matrix}$ according to formula (14) to formula (16), since the quadratic term coefficients of the second derivative and the first derivative are both larger than 0, as |İ_(n) ^(θ) ^(n) ^(*)| changes from 0 to |Δİ_(n) ^(θ) ^(n) ^(*,max)|, ΔI_(Σ) has at most one extreme point which is a minimum point; ΔI_(Σ) is maximum when Δİ_(n) ^(θ) ^(n) ^(*)=0 or Δİ_(n) ^(θ) ^(n) ^(*)=Δİ_(n) ^(θ) ^(n) ^(*,max); at least one Δİ_(ξ) ^(n) inevitably exists in the sector shadow area so that: |Δİ* ₁ +Δİ* ₂ + . . . +Δİ* _(n-1) +Δİ _(n) ^(ξ) |>|Δİ* ₁ +Δİ* ₂ + . . . +Δİ* _(n-1)+0|  (17) when the overall error of superimposed multi-terminal current errors is maximum, Δİ*_(n)=Δİ_(n) ^(θ) ^(n) ^(*,max); numbering four boundaries of the sector shadow area corresponding to Δİ_(n), the boundaries 1 and 2 are straight boundaries, and the boundaries 3 and 4 are arc boundaries; wherein Δİ_(I) ending on the boundary 1 can be decomposed into: Δİ _(n) =Δİ* _(n,x) +Δİ _(n,y)  (18) wherein Δİ*_(n,x) is the projection of Δİ_(n) on the vertical line of the boundary 1, and Δİ_(n,y) is the projection of Δİ_(n) on the boundary 1; when Δİ*_(n,y) ends near the A₁ side, ΔI_(Σ.max) ² is: $\begin{matrix} \begin{matrix} {{\Delta I_{\Sigma.\max}^{2}} = {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n}^{*}}}❘}^{2}} \\ {= {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}} + {\Delta{\overset{.}{I}}_{n,y}^{*}}}❘}^{2}} \\ {= {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}^{2} + {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘}^{2} - {2 \cdot {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘} \cdot}}} \\ {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}{\cos\left( {\frac{\pi}{2} - \delta_{n}^{*}} \right)}} \end{matrix} & (19) \end{matrix}$ wherein δ*_(n) is the phase difference between Δİ*₁+Δİ*₂+ . . . +Δİ*_(n-1) and Δİ*_(n,x), and δ*_(n)∈[0, 2π); since Δİ*₁, Δİ*₂, . . . , Δİ*_(n-1), Δİ*_(n,x), Δİ_(n,y) are independent of each other, the derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is: $\begin{matrix} {\frac{d\left( {\Delta I_{\Sigma}^{2}} \right)}{d\left( {❘{\Delta{\overset{.}{I}}_{n,y}}❘} \right)} = {{2{❘{\Delta{\overset{.}{I}}_{n,y}}❘}} - {2{{\cos\left( {\frac{\pi}{2} - \delta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}}}} & (20) \end{matrix}$ the second derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is: $\begin{matrix} {\frac{d^{2}\left( {\Delta I_{\Sigma}^{2}} \right)}{{d\left( {❘{\Delta{\overset{.}{I}}_{n,y}}❘} \right)}^{2}} = {2 > 0}} & (21) \end{matrix}$ when Δİ*_(n,y) ends near the A₂ side, ΔI_(Σ.max) ² is: $\begin{matrix} \begin{matrix} {{\Delta I_{\Sigma.\max}^{2}} = {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n}^{*}}}❘}} \\ {= {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}} + {\Delta{\overset{.}{I}}_{n,y}^{*}}}❘}^{2}} \\ {= {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}^{2} + {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘}^{2} - {2 \cdot {❘{\Delta{\overset{.}{I}}_{n,y}^{*}}❘} \cdot}}} \\ {{❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}{\cos\left( {\frac{3\pi}{2} - \delta_{n}^{*}} \right)}} \end{matrix} & (22) \end{matrix}$ at this moment, the derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is: $\begin{matrix} {\frac{d\left( {\Delta I_{\Sigma}^{2}} \right)}{d\left( {❘{\Delta{\overset{.}{I}}_{n,y}}❘} \right)} = {{2{❘{\Delta{\overset{.}{I}}_{n,y}}❘}} - {2{{\cos\left( {\frac{3\pi}{2} - \delta_{n}^{*}} \right)} \cdot {❘{{\Delta{\overset{.}{I}}_{1}^{*}} + {\Delta{\overset{.}{I}}_{2}^{*}} + \ldots + {\Delta{\overset{.}{I}}_{n - 1}^{*}} + {\Delta{\overset{.}{I}}_{n,x}^{*}}}❘}}}}} & (23) \end{matrix}$ when Δİ*_(n,y) ends near the A₂ side, the second derivative of ΔI_(Σ) ² with respect to |Δİ_(n,y)| is the same as formula (21): denoting Δİ_(n,y) ending at A₁ as Δİ_(n,y) ^(A1), and denoting Δİ_(n,y) ending at A₂ as Δİ_(n,y) ^(A2); respectively denoting derivatives in formula (20) and formula (23) as f_(n,y) ¹ and f_(n,y) ².
 5. The multi-terminal differential protection setting method for a distributed generator T-connected distribution network according to claim 1, wherein step S7 of solving the optimization problem of the non-convex constraint space comprises the following specific steps: step S701: inputting data matrixes of İ_(1.tr)-İ_(n.tr), λ₁-λ_(n) and α₁-α_(n); $\begin{matrix} \begin{bmatrix} I_{1.{tr}} & \beta_{1,{tr}} & \lambda_{1} & \alpha_{1} \\ I_{2.{tr}} & \beta_{2,{tr}} & \lambda_{2} & \alpha_{2} \\  \vdots & \vdots & \vdots & \vdots \\ I_{n,{tr}} & \beta_{n,{tr}} & \lambda_{n} & \alpha_{n} \end{bmatrix} & (28) \end{matrix}$ wherein I_(1.tr)-I_(n.tr) are amplitudes of İ_(1.tr)-İ_(n.tr), and β_(1.tr)-β_(1.tr) are phases of İ_(1.tr)-İ_(n.tr); wherein the specific values of λ₁-λ_(n) and α₁-α_(n) can be determined by the type of a current transformer and the data synchronization accuracy of a data synchronization method, and İ_(1.tr)-İ_(n.tr) are fault current flowing through primary and secondary fusing intelligent switches when a fault outside the protected area is most severe; step S702: determining to discretize the constraint space of the optimization problem into (4+2m)n points according to λ₁λ_(n) and α₁-α_(n); step S703: generating a discrete point matrix; $\begin{matrix} \begin{bmatrix} {{I_{1.{tr}}\left( {1 - \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \alpha_{1}} \right)}} & {I_{2.{tr}}\left( {1 - \lambda_{2}} \right)\angle\left( {\beta_{2.{tr}} - \alpha_{2}} \right)} & \ldots & {I_{n.{tr}}\left( {1 - \lambda_{n}} \right)\angle\left( {\beta_{n.{tr}} - \alpha_{n}} \right)} \\ {I_{1.{tr}}\left( {1 - \lambda_{1}} \right)\angle\left( {\beta_{1.{tr}} + \alpha_{1}} \right)} & {I_{2.{tr}}\left( {1 - \lambda_{2}} \right)\angle\left( {\beta_{2.{tr}} + \alpha_{2}} \right)} & \ldots & {I_{n.{tr}}\left( {1 - \lambda_{n}} \right)\angle\left( {\beta_{n.{tr}} + \alpha_{n}} \right)} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \alpha_{1}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \alpha_{2}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \alpha_{n}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \alpha_{1}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \alpha_{2}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \alpha_{n}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \frac{\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \frac{\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \frac{\alpha_{n}}{2^{m}}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \frac{3\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \frac{3\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \frac{3\alpha_{n}}{2^{m}}} \right)}} \\  \vdots & \vdots & \ldots & \vdots \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} - \frac{\left( {{2k} + 1} \right)\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} - \frac{\left( {{2k} + 1} \right)\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} - \frac{\left( {{2k} + 1} \right)\alpha_{n}}{2^{m}}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \frac{\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \frac{\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \frac{\alpha_{n}}{2^{m}}} \right)}} \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \frac{3\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \frac{3\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \frac{3\alpha_{n}}{2^{m}}} \right)}} \\  \vdots & \vdots & \ldots & \vdots \\ {{I_{1.{tr}}\left( {1 + \lambda_{1}} \right)}{\angle\left( {\beta_{1.{tr}} + \frac{\left( {{2k} + 1} \right)\alpha_{1}}{2^{m}}} \right)}} & {{I_{2.{tr}}\left( {1 + \lambda_{2}} \right)}{\angle\left( {\beta_{2.{tr}} + \frac{\left( {{2k} + 1} \right)\alpha_{2}}{2^{m}}} \right)}} & \ldots & {{I_{n.{tr}}\left( {1 + \lambda_{n}} \right)}{\angle\left( {\beta_{n.{tr}} + \frac{\left( {{2k} + 1} \right)\alpha_{n}}{2^{m}}} \right)}} \end{bmatrix} & (29) \end{matrix}$ wherein k≥1 and 2k−1≤2^(m)−1 step S704: for 1-n columns in the matrix in step S703, performing the following (4+2^(m))^(n) operations respectively for different points x_({dot over (ι)}1,1), x_({dot over (ι)}2,2), . . . , x_({dot over (ι)}n,n) of 4+2^(m) rows; ΔI _(Σ) =|x _({dot over (ι)}1,1) +x _({dot over (ι)}2,2) + . . . +x _({dot over (ι)}n,n)|  (30) wherein 1≤{dot over (ι)}1, {dot over (ι)}2, . . . , {dot over (ι)}_(n)≤4+2^(m); step S705: searching for the maximum value in the (4+2^(m))^(n) results calculated in step S704, i.e., the approximate solution ΔI_(Σ.max)′ of a global optimal solution, and approximately replacing ΔI_(Σ.max) by ΔI_(Σ.max)′ as the protection setting value in formula (9).
 6. The multi-terminal differential protection setting method for a distributed generator T-connected distribution network according to claim 1, wherein step S8 of analyzing the solution result of the optimization problem of the non-convex constraint space comprises: setting up a simulation model in PSCAD/EMTDC simulation software by using the system topology of a high-proportion distributed generator T-connected distribution network with a power electronic transformer as a transmission and distribution interface.
 7. The multi-terminal differential protection setting method for a distributed generator T-connected distribution network according to claim 1, wherein step S8 of verifying the solution result of the optimization problem of the non-convex constraint space comprises: conducting contrast verification through results of sector boundary traversal and six-point, eight-point and twelve-point fast traversal methods. 